# The script plots the grid at year 20 where color is linked to genetic distance
# uses PCA to shrink the 244 dimensions of microsat loci into 3 dimensions of red, green, blue.

# Plot the hexagonal grid at year 20
library(plotrix)

for(i in 1:1){
	# Read the patterns
	patterns <- read.csv(paste("sweep",i,"/sweep",i,".grid.patterns",sep=""),header=FALSE)
	# PCA them
	pat <- patterns[,4:247]
	# PCA
	pat_pca <- prcomp(pat)
	red_col <- predict(pat_pca)[,1]
	grn_col <- predict(pat_pca)[,2]
	blu_col <- predict(pat_pca)[,3]
	red_col <- red_col+abs(min(red_col))
	red_col <- red_col/max(red_col)
	grn_col <- grn_col+abs(min(grn_col))
	grn_col <- grn_col/max(grn_col)
	blu_col <- blu_col+abs(min(blu_col))
	blu_col <- blu_col/max(blu_col)
        crypts <- read.csv(paste("sweep",i,"/sweep",i,".crypts",sep=""))
	# A clone is a unique ms pattern
	# crypt_id, pattern_id, cid, year, gen
        clones <- read.csv(paste("sweep",i,"/sweep",i,".grid",sep=""))

	clonecolors <- NULL
	clonecount <- max(clones[,2])
	# Ordered pattern_id 1..# of patterns assigned color
	clonecolors <- cbind(c(1:clonecount),red_col,grn_col,blu_col)

for(j in 1:39){
	# take time slice
	clonestimepoint <- clones[clones[,4]==j*182,]
	# plot the matrix
	mat <- matrix(crypts[,1],nrow=256,ncol=256,byrow=TRUE)
	matcol <- mat
	for(x in 1:256){
	for(y in 1:256){
	if(x==256&y==256){
		matcol[x,y]<-rgb(255,255,255,max=255)
	} else {
matcol[x,y]<-rgb(clonecolors[clonestimepoint[(x-1)*256+y,2],2],clonecolors[clonestimepoint[(x-1)*256+y,2],3],clonecolors[clonestimepoint[(x-1)*256+y,2],4])
	}
	}
	}
png(paste("sweep",i,"/sweep",i,".time.",j,".grid.png",sep=""),width=512,height=512,units="px",res=72)
color2D.matplot(mat,cellcolors=matcol,xlab=NA,ylab=NA,do.hex=TRUE,border=NA,axes=FALSE)
dev.off()

	print (paste("Done sweep ",j*182))
}

for(j in 41:41){
	# take time slice
	clonestimepoint <- clones[clones[,4]==7300,]
	# plot the matrix
	mat <- matrix(crypts[,1],nrow=256,ncol=256,byrow=TRUE)
	matcol <- mat
	for(x in 1:256){
	for(y in 1:256){
	if(x==256&y==256){
		matcol[x,y]<-rgb(255,255,255,max=255)
	} else {
matcol[x,y]<-rgb(clonecolors[clonestimepoint[(x-1)*256+y,2],2],clonecolors[clonestimepoint[(x-1)*256+y,2],3],clonecolors[clonestimepoint[(x-1)*256+y,2],4])
	}
	}
	}
png(paste("sweep",i,"/sweep",i,".time.",j,".grid.png",sep=""),width=512,height=512,units="px",res=72)
color2D.matplot(mat,cellcolors=matcol,xlab=NA,ylab=NA,do.hex=TRUE,border=NA,axes=FALSE)
dev.off()

	print (paste("Done sweep ",7300))
}




}


png(paste("sweep",i,"/sweep",i,".grid.png",sep=""),width=2.5,height=2.5,units="in",res=72)
color2D.matplot(mat,cellcolors=matcol,xlab=NA,ylab=NA,do.hex=TRUE,border=NA,axes=FALSE,mfrow=c(1,1),pty="s",mar=c(0,0,0,0),mai=c(0,0,0,0))
dev.off()

